Raport dot. elektrod grafenowych

1 Podsumowanie

Niniejszy raport przedstawia szczegółówą analizę właściwości elektrod grafenowych na podstawie zbioru danych. Dane zawierają 925 próbek, z czego znaczna jej część jest niekompletna. Do ich uzupełnienia zastosowano metodę k-nearest neighbors. Następnie przeprowadzono analizę statystyczną, sprawdzono rozkłady wartości atrybutów oraz, na podstawie macierzy korelacji, znaleziono kilka ciekawych zależności.

Jako rezultat przeprowadzonej analizy zaobserwowano:
1. Udział węgla ma istotny wpływ na właściwości materiału. Zwiększenie zawartości węgla w materiale wysoce wpływa na zmniejszoną pojemność materiału. W zamian zyskujemy zwiększenie okna stabilności.
2. Pojemność rośnie wraz z udziałem tlenu. Niski udział tlenu (<25%) praktycznie gwarantuje niską pojemność materiału, ale zapewnia wysokie okno stabilności. Udziął azotu ma natomiast marginalny wpływ na pojemność.
3. Pojemność jest negatywnie skorelowana z oknem stabilności. Trudno jest znaleźć materiał który ma jednocześnie wysoką stabilność i ma dużą pojemność.

W części modelowania podjęto próbę przewidzenia pojemności materiału na podstawie pozostałych informacji ze zbioru. Skorzystano z dwóch modeli uczenia maszynowego:
1. Model regresji liniowej, który okazał się niewystarczający, szczególnie dla próbek o wysokiej pojemności,
2. Model Random Forest, który poradził sobie zdecydowanie lepiej, dla niższych wartości pojemności trafia praktycznie w punkt, a dla wyższych tylko nieznacznie je niedoszacowuje.

Na zakończenie zastosowano analizę SHAP dla modelu Random Forest pokazując przykładowe rozumowanie algorytmu.

2 Przygotowanie zbioru do analizy

2.1 Wykorzystane biblioteki

Do przygotowania raportu skorzystano z następujących bibliotek przedstawionych na zajęciach:

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggsci)
library(VIM)
library(tibble)
library(knitr)
library(plotly)
library(randomForest)
library(caret)
library(shapper)
library("DALEX2")

2.2 Zapewnienie powtarzalności wyników.

Dla zapewnienia powtarzalności wyników skorzystano z funkcji set.seed()

set.seed(123)

2.3 Wczytanie danych

Zbiór danych dotyczy właściwości elektrod grafenowych i zawiera szczegółowe informacje na temat ich budowy, warunków testowania oraz uzyskanej wydajności.

Plik jest wczytywany przy założeniu, że jest umieszczony w tym samym katalogu co plik z tym raportem.

df <- read.csv('data.csv')
df <- as_tibble(df)

Wszystkie kolumny zostały przemianowane na jej skrócone odpowiedniki, aby ułatwić przetwarzanie danych. Usunięto kolumnę ref, ponieważ jest to jedynie identyfikator którego nie będziemy używać przy analizie danych. Dodatkowo, usuwamy pw_limits, z uwagi na to, że powiela ona informację z pw_low i pw_up

df <- df %>% 
  rename(
    ref = `Ref.`,
    pw_limits = `Limits.of.Potential.Window..V.`,
    pw_low = `Lower.Limit.of.Potential.Window..V.`,
    pw_up = `Upper.Limit.of.Potential.Window..V.`,
    pw = `Potential.Window..V.`,
    current_density = `Current.Density..A.g.`,
    capacitance = `Capacitance..F.g.`,
    ssa = `Specific.Surface.Area..m.2.g.`,
    rct = `Charge.Transfer.Resistance..Rct...ohm.`,
    rs = `Equivalent.Series.Resistance..Rs...ohm.`,
    electrode_conf = `Electrode.Configuration`,
    pore_size = `Pore.Size..nm.`,
    pore_volume = `Pore.Volume..cm.3.g.`,
    id_ig_ratio = `Ratio.of.ID.IG`,
    n_at = `N.at.`,
    c_at = `C.at.`,
    o_at = `O.at.`,
    electrolyte_formula = `Electrolyte.Chemical.Formula`,
    electrolyte_cond = `Electrolyte.Ionic.Conductivity`,
    electrolyte_conc = `Electrolyte.Concentration..M.`,
    cell_conf = `Cell.Configuration..three.two.electrode.system.`
  ) %>%
  select(-ref, -pw_limits)
head(df)

## # A tibble: 6 × 19
##   pw_low pw_up    pw current_density capacitance   ssa   rct    rs
##    <dbl> <dbl> <dbl>           <dbl>       <dbl> <dbl> <dbl> <dbl>
## 1      0   0.8   0.8               1         680  186.  NA    7.7 
## 2      0   1     1                 1         367  537    6.1  1.95
## 3      0   1     1                 2         338  537    6.1  1.95
## 4      0   1     1                 5         283  537    6.1  1.95
## 5      0   1     1                10         246  537    6.1  1.95
## 6      0   0.5   0.5               1         872  168.  NA    0.8 
## # ℹ 11 more variables: electrode_conf <chr>, pore_size <dbl>,
## #   pore_volume <dbl>, id_ig_ratio <dbl>, n_at <dbl>, c_at <dbl>, o_at <dbl>,
## #   electrolyte_formula <chr>, electrolyte_cond <int>, electrolyte_conc <dbl>,
## #   cell_conf <chr>

2.4 Przetwarzanie brakujących danych

Poniżej przedstawiona została ilość brakujących danych dla każdej kolumny.

count(df)

## # A tibble: 1 × 1
##       n
##   <int>
## 1   925

df %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  as.data.frame()

##   pw_low pw_up pw current_density capacitance ssa rct  rs electrode_conf
## 1      4     4  5              16          17 572 786 772              0
##   pore_size pore_volume id_ig_ratio n_at c_at o_at electrolyte_formula
## 1       769         729         596  690  699  703                   0
##   electrolyte_cond electrolyte_conc cell_conf
## 1               99               62         0

Dodatkowo, później przy analizie danych zauważono, że niektóre wartości kolumn tekstowych są puste, lecz w formie pustej treści a nie jako “NA”. Takie przypadki również podmienimy na NA.

df <- df %>%
    mutate(across(
      .col = where(is.character),
      .fns = ~na_if(.x, '')
    ))

Skala brakujących danych jest poważna. W kilku kolumnach zdecydowana większość obserwacji ma wartości puste. Jako rozwiązanie zaproponowano zastosowanie algorytmu K-nearest-neighbors, aby uzupełnić kolumny z jak najmniejszym możliwym wpływem na statystyki kolumny jak np. mediana wartości.

Dodatkowo, dla kilku wybranych metryk dodamy nową kolumnę zapisującą informację czy dana wartość była początkowo pusta. Zrobiono to, ponieważ bywają sytuację w których brak informacji również jest wartościową informacją (przykładem - stosunek odmów odpowiedzi w ankietach wyborczych jest różny dla sympatyków różnych partii).

df <- df %>%
  mutate(rct_na = is.na(rct), rs_na = is.na(rs), n_at_na = is.na(n_at))

df <- kNN(df, k=5, imp_var = FALSE)

3 Analiza zbioru

3.1 Informacje o zbiorze

Dane po przetworzeniu brakujących danych zawiewrają następujące dane:

• pw_low, pw_up (V): Przedział, w którym testowany materiał pozostawał stabilny.
• current_density (A/g): Ilość napięcia na jednostkę masy.
• capacitance (F/g): Pojemność właściwa, czyli kluczowy parametr wyjściowy określający zdolność materiału elektrody do magazynowania ładunku na jednostkę masy (w faradach na gram).
• ssa (m²/g): Powierzchnia właściwa (SSA), czyli całkowita powierzchnia materiału elektrody na jednostkę masy (w metrach kwadratowych na gram).
• rct (Rct, Ω): Opór transferu ładunku pomiędzy elektrodą a elektrolitem.
• rs (Rs, Ω): Całkowity opór wewnętrzny superkondensatora.
• electrode_conf: Konfiguracja elektrody, opisująca skład chemiczny i strukturę materiału kompozytowego użytego w elektrodzie.
• pore_size (nm): Rozmiar porów w materiale elektrody (w nanometrach).
• pore_volume (cm³/g): Całkowita objętość porów na jednostkę masy materiału (w centymetrach sześciennych na gram).
• id_ig_ratio: Wskaźnik poziomu defektów i nieuporządkowania w strukturze grafenu.
• n_at (%): Procent atomów azotu.
  
• c_at (%): Procent atomów węgla.
  
• o_at (%): Procent atomów tlenu.
  
• electrolyte_formula: Wzór chemiczny elektrolitu użytego w badaniach.
• electrolyte_cond: Wartość wskazująca jak dobrze elektrolit przewodzi prąd.
• electrolyte_conc (M): Stężenie molowe elektrolitu.
• cell_conf: Konfiguracja ogniwa pomiarowego, określająca, czy w badaniach zastosowano układ dwu- czy trójelektrodowy.
• rct_na rs_na n_at_na: Wartość boolean, opisująca czy kolumny rct, rs, n_at były początkowo puste.

Zbiór zawiera 925 próbek.

head(df)

##   pw_low pw_up  pw current_density capacitance   ssa  rct   rs
## 1      0   0.8 0.8               1         680 186.3 18.5 7.70
## 2      0   1.0 1.0               1         367 537.0  6.1 1.95
## 3      0   1.0 1.0               2         338 537.0  6.1 1.95
## 4      0   1.0 1.0               5         283 537.0  6.1 1.95
## 5      0   1.0 1.0              10         246 537.0  6.1 1.95
## 6      0   0.5 0.5               1         872 168.2  0.6 0.80
##                     electrode_conf pore_size pore_volume id_ig_ratio n_at  c_at
## 1                   CNF/RGO/moOxNy      3.45       0.046        1.45 2.10 18.45
## 2 sulfur-doped graphene foam (SGF)     37.00       0.500        1.28 0.00 85.60
## 3 sulfur-doped graphene foam (SGF)     37.00       0.500        1.28 0.00 85.60
## 4 sulfur-doped graphene foam (SGF)     37.00       0.500        1.28 0.00 85.60
## 5 sulfur-doped graphene foam (SGF)     37.00       0.500        1.28 0.00 85.60
## 6                      moS2-moO2/G      3.45       0.170        1.22 1.88 30.75
##    o_at electrolyte_formula electrolyte_cond electrolyte_conc
## 1 44.31               H2SO4                7                1
## 2  9.10                 KOH                6                6
## 3  9.10                 KOH                6                6
## 4  9.10                 KOH                6                6
## 5  9.10                 KOH                6                6
## 6 52.48                 KOH                6                2
##                cell_conf rct_na rs_na n_at_na
## 1 three-electrode system   TRUE FALSE   FALSE
## 2   two-electrode system  FALSE FALSE   FALSE
## 3   two-electrode system  FALSE FALSE   FALSE
## 4   two-electrode system  FALSE FALSE   FALSE
## 5   two-electrode system  FALSE FALSE   FALSE
## 6 three-electrode system   TRUE FALSE    TRUE

nrow(df)

## [1] 925

Poniżej wyświetlono kilka metryk dla każdej kolumny ze zbioru

summary(df)

##      pw_low            pw_up               pw         current_density 
##  Min.   :-1.1000   Min.   :-0.2000   Min.   :0.4000   Min.   :  0.05  
##  1st Qu.:-0.3000   1st Qu.: 0.4000   1st Qu.:0.6000   1st Qu.:  1.00  
##  Median : 0.0000   Median : 0.6000   Median :0.8000   Median :  2.00  
##  Mean   :-0.2335   Mean   : 0.6302   Mean   :0.8628   Mean   :  5.80  
##  3rd Qu.: 0.0000   3rd Qu.: 0.8000   3rd Qu.:1.0000   3rd Qu.:  5.00  
##  Max.   : 0.2000   Max.   : 3.5000   Max.   :3.5000   Max.   :200.00  
##   capacitance          ssa                rct               rs        
##  Min.   :   1.4   Min.   :   8.896   Min.   : 0.080   Min.   : 0.200  
##  1st Qu.: 150.0   1st Qu.: 105.000   1st Qu.: 1.100   1st Qu.: 0.500  
##  Median : 262.5   Median : 159.970   Median : 2.010   Median : 0.760  
##  Mean   : 419.3   Mean   : 427.082   Mean   : 2.761   Mean   : 1.351  
##  3rd Qu.: 509.9   3rd Qu.: 577.100   3rd Qu.: 3.520   3rd Qu.: 2.000  
##  Max.   :3344.1   Max.   :2400.000   Max.   :24.200   Max.   :17.500  
##  electrode_conf       pore_size      pore_volume      id_ig_ratio  
##  Length:925         Min.   : 0.53   Min.   :0.0200   Min.   :0.12  
##  Class :character   1st Qu.: 6.20   1st Qu.:0.0460   1st Qu.:0.94  
##  Mode  :character   Median : 6.20   Median :0.2160   Median :1.09  
##                     Mean   :10.29   Mean   :0.3702   Mean   :1.09  
##                     3rd Qu.:10.80   3rd Qu.:0.4610   3rd Qu.:1.18  
##                     Max.   :44.13   Max.   :2.3500   Max.   :2.90  
##       n_at            c_at            o_at       electrolyte_formula
##  Min.   : 0.00   Min.   : 1.40   Min.   : 1.90   Length:925         
##  1st Qu.: 1.88   1st Qu.:30.75   1st Qu.:12.83   Class :character   
##  Median : 1.88   Median :54.18   Median :18.50   Mode  :character   
##  Mean   : 2.49   Mean   :55.97   Mean   :22.49                      
##  3rd Qu.: 2.10   3rd Qu.:81.30   3rd Qu.:26.70                      
##  Max.   :23.82   Max.   :98.10   Max.   :54.28                      
##  electrolyte_cond electrolyte_conc  cell_conf           rct_na       
##  Min.   :1.000    Min.   :0.100    Length:925         Mode :logical  
##  1st Qu.:6.000    1st Qu.:1.000    Class :character   FALSE:139      
##  Median :6.000    Median :1.000    Mode  :character   TRUE :786      
##  Mean   :5.818    Mean   :2.698                                      
##  3rd Qu.:7.000    3rd Qu.:6.000                                      
##  Max.   :8.000    Max.   :6.000                                      
##    rs_na          n_at_na       
##  Mode :logical   Mode :logical  
##  FALSE:153       FALSE:235      
##  TRUE :772       TRUE :690      
##                                 
##                                 
## 

3.2 Analiza wartości atrybutów

Poniżej umieszczono rozkłady wartości wszystkich wartości numerycznych

Widzimy, że wartości w zbiorze są szeroko rozproszone, a często pojawia się w nich szum informacyjny, gdzie maksymalnie kilka wartości odbiega od ogólnego trendu. Spróbujmy odfiltrować ten szum o 10% najniższych i najwyższych wartości.

Oprócz kolumn z wartościami numerycznymi, możemy też skategoryzować atrybut “cell_conf” i policzyć jego występowanie. Podobnie można zrobić dla “electrode_conf” i “electrolyte_formula”, lecz liczba ich kategorii jest zbyt duża aby miało sens wyświetlenie ich na wykresach.

3.3 Analiza korelacji między zmiennymi

Do znalezienia korelacji między zmiennymi wykorzystano macierz korelacji.

Wyszukajmy wszystkie koreleacje najbardziej warte odnotowania. Jako takie uznamy takie pary, których korelacja znajduje się między 0.99 i 0.4 oraz między -0.99 i -0.4.

cor_df_high <- cor_df %>%
  filter((abs(correlation) > 0.4 & abs(correlation) < 0.99)) %>%
  filter(var1 > var2) %>%
  arrange(correlation)
kable(cor_df_high)

var1	var2	correlation
o_at	c_at	-0.7010226
capacitance	c_at	-0.4553412
pw	o_at	-0.4164487
pw	c_at	0.4714214
pw_up	pw	0.5970643
pw_up	pw_low	0.6370278

Analizując znalezione korelacje możemy pominąć korelacje pw_up z pw_low oraz pw_up z pw, gdyż wydają się dosyć naturalne. Jako najciekawsze uznałem zależności między: - capacitance a c_at, - pw a c_at.

3.4 Analiza trendów

Obydwa wykresy pokazują ciekawe właściwości węgla. Możemy z nich wyczytać, że wraz ze zwiększeniem zawartości węgla w materiale jego pojemność właściwa spada. Z drugiej strony rośnie okno w którym materiał pozostaje stabilny. Sprawdźmy jeszcze zależność między pojemnością a oknem stabilności, aby sprawdzić czy one też nie są ze sobą skorelowane.

Jak widać na wykresie, korelacja między nimi również występuje. Natomiast nakreślony trend dla powyżej 2000 F/g jest bardzo wątpliwy z uwagi na stałą wartość 0.5V dla tych pojemności, która może wynikać z metody uzupełniania pustych wartości lub samego zbioru danych.

Przyjrzyjmy się jeszcze pojemności w relacji do składu atomów.

Z wykresu wyżej widzimy, że większy udział tlenu wspiera do pewnego stopnia pojemność materiału. Tak jak wcześniej zauważyliśmy, wpływ węgla jest negatywny na pojemność, a jego odbicie na końcu prawdopodobnie jest jedynie błędem z powodu małej próbki. Ostatnim wnioskiem jest bardzo niski udział azotu we wszystkich materiałach, a jego zmiana ma marginalny wpływ na pojemność.

3.5 Model przewidywania cech materiałów

W tej części spróbujemy przewidzieć pojemność materiału na podstawie tego zbioru danych.Na początku spróbujmy znaleźć wartościowe parametry. Skorzystamy z poprzedniego kodu. Przyjmijmy, że kolumny warte zachowania to są te których korelacja jest powyżej 0.1 lub poniżej -0.1 (z pominięciem pw_low i pw_up, gdyż zakładam ich dużą korelację z pw).

Spróbujmy zamodelować pojemność regresją liniową.

cor_df_cap <- cor_df %>%
  filter(var1 == 'capacitance' & var2 != 'capacitance' & abs(correlation) > 0.1) %>%
  arrange(correlation)
kable(cor_df_cap)

var1	var2	correlation
capacitance	c_at	-0.4553412
capacitance	pw	-0.3384832
capacitance	ssa	-0.2209595
capacitance	pw_up	-0.1047333
capacitance	pore_volume	-0.1010717
capacitance	id_ig_ratio	0.1026857
capacitance	pw_low	0.1752674
capacitance	o_at	0.3204094

df_capacitance <- df %>%
  select(capacitance, c_at, pw, ssa, pore_volume, id_ig_ratio, o_at)

inTraining <- 
    createDataPartition(
        y = df_capacitance$capacitance,
        p = .75,
        list = FALSE)

training <- df_capacitance[ inTraining,]
testing  <- df_capacitance[-inTraining,]

ctrl <- trainControl(
    method = "repeatedcv",
    number = 2,
    repeats = 5
)

fit <- train(capacitance ~ .,
             data = training,
             method = "lm",
             trControl = ctrl
)
rfClasses <- predict(fit, newdata = testing)

metrics <- postResample(pred = rfClasses, obs = testing$capacitance)
metrics

##        RMSE    Rsquared         MAE 
## 387.8635335   0.2515056 268.3858543

results_df <- data.frame(
  Actual = testing$capacitance,
  Predicted = rfClasses
)

min_val <- min(c(results_df$Actual, results_df$Predicted))
max_val <- max(c(results_df$Actual, results_df$Predicted))

p <- ggplot(results_df, aes(x = Actual, y = Predicted)) +
  geom_point(alpha = 0.6, color = "#1f77b4") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red", linewidth = 1) +
  geom_smooth(method = "lm", color = "darkgreen", se = FALSE) +
  coord_fixed(xlim = c(min_val, max_val), ylim = c(min_val, max_val)) +
  labs(
    title = "Prognozowana vs. rzeczywista pojemność - regresja liniowa",
    x = "Rzeczywista pojemność (F/g)",
    y = "Przewidywana pojemność (F/g)"
  ) +
  theme_minimal()

ggplotly(p)

Na wykresie widzimy 2 linie - zielona to nasz wytrenowany model, a czerwona to idealna regresja dla tego zbioru danych. Jak można na wykresie zauważyć, nasz model kompletnie nie radzi sobie z większymi pojemnościami. Spróbujmy skorzystać z algorytmu Random Forest Regression.

head(training)

##    capacitance  c_at  pw   ssa pore_volume id_ig_ratio  o_at
## 2          367 85.60 1.0 537.0        0.50        1.28  9.10
## 3          338 85.60 1.0 537.0        0.50        1.28  9.10
## 4          283 85.60 1.0 537.0        0.50        1.28  9.10
## 6          872 30.75 0.5 168.2        0.17        1.22 52.48
## 7          143 54.18 0.6 385.0        0.10        1.00 12.83
## 10         483 54.18 0.6 385.0        0.31        1.21 12.83

rf_fit <- train(capacitance ~ .,
               data = training,
               method = "rf",
               trControl = ctrl,
               tuneLength = 3
)
rf_predicts <- predict(rf_fit, newdata = testing)
rf_metrics <- postResample(pred = rf_predicts, obs = testing$capacitance)

results_df <- data.frame(
  Actual = testing$capacitance,
  Predicted = rf_predicts
)

min_val <- min(c(results_df$Actual, results_df$Predicted))
max_val <- max(c(results_df$Actual, results_df$Predicted))

p <- ggplot(results_df, aes(x = Actual, y = Predicted)) +
  geom_point(alpha = 0.6, color = "#1f77b4") +
  geom_abline(intercept = 0, slope = 1, linetype = "solid", color = "red", linewidth = 1) +
  geom_smooth(method = "lm", color = "darkgreen", se = FALSE, linewidth = 0.8) +
  coord_fixed(xlim = c(min_val, max_val), ylim = c(min_val, max_val)) +
  labs(
    title = "Prognozowana vs. rzeczywista pojemność - RF",
    x = "Rzeczywista pojemność (F/g)",
    y = "Przewidywana pojemność (F/g)"
  ) +
  theme_minimal()

ggplotly(p)

Jak widzimy na wykresie, algorytm Random Forest poradził sobie o rząd wielkości lepiej aniżeli regresja liniowa. Spójrzmy na sam koniec jak wyglądają przykładowe wartości Shapleya dla jednego przykładu.

testing[1,]

##   capacitance  c_at  pw   ssa pore_volume id_ig_ratio  o_at
## 1         680 18.45 0.8 186.3       0.046        1.45 44.31

ive_rf <- individual_variable_effect(rf_fit, data = testing,
                                     new_observation = testing[1,])
plot(ive_rf)